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Abstract 

We describe a Godunov-type magnetohydrodynamic (MHD) code based on the Miyoshi & Kusano (2005) 
solver which can be used to solve various astrophysical hydrodynamic and MHD problems. The energy 
equation is in the form of entropy conservation. The code has been implemented on several different 
coordinate systems: 2.5D axisymmetric cylindrical coordinates, 2D Cartesian coordinates, 2D plane polar 
coordinates, and fully 3D cylindrical coordinates. Viscosity and diffusivity are implemented in the code to 
control the accretion rate in the disk and the rate of penetration of the disk matter through the magnetic 
field lines. The code has been utilized for the numerical investigations of a number of different astrophysical 
problems, several examples of which are shown. 
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1. Introduction 


A number of numerical magnetohydrodynamic 
(MHD) codes have been developed for model¬ 
ing plasma flows in astrophysics. Some of the 


man 

1992), FLASH (Fryxell et al. 

2000), PLUTO 

(Mignone et al. 

2007), and ATHENA (Stone et 

al. 

00 
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Skinner & OstrikerJ 2010). A num- 


her of different numerical algorithms have also 
been developed for the numerical integration of 
the MHD equations including different approaches 


for the spatial and temporal approximations (Brio 

& Wu, 1988; Cockburn, Lin & Shu 

1989; 

Dai & 

Woodward 1994a b 

Ryu & Jones 1995; Balsara & 

Spicer 

1999; Gurski 

2004; Ustyugov et al. 

2009), 


ims for the approximate solu- 


and different algorit 
tion of the Riemann problem (Brio & Wu[ 1988 


Li| [2005j |Miyoshi & Kusano| |2005[ |Miyoshi et al. 

2010D. 


In this work, we describe a code developed by 
our group for the numerical modeling of astrophys¬ 


ical MHD flows. The code has been developed for 
use in several coordinate systems: (1) 2.5D axisym¬ 
metric cylindrical coordinates (r,z); (2) 2D polar 
coordinates (r, <f)] (3) 3D cylindrical coordinates 
( r,(j>,z ); and (4) a Cartesian (x, v) geometry which 
is used to conduct tests of the ideal and non-ideal 
MHD modules. The difference between our code 
and the above-mentioned codes lies in the specifics 
of the astrophysical problems that we solve. Firstly, 
in the astrophysical regimes that we study, strong 
shocks (where the energy dissipation cannot be ne¬ 
glected) are not expected to occur. This permits the 
use of the entropy conservation equation instead 
of the full energy equation. The advantage of this 
approach is that the entropy conservation equation 
does not contain terms that differ significantly in 
magnitude. For example, in the energy equation, 
the largest terms are the gravitational energy and 
the kinetic energy of the azimuthal motion; these 
can be much larger than the internal energy and 
the energy of the poloidal motion near a gravitat¬ 
ing body (for example, a planet or star). Addition- 
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ally, in the vicinity of a magnetized star, the mag¬ 
netic energy density can significantly exceed the 
matter pressure (and the energy-density of thermal 
energy), leading to significant errors when comput¬ 
ing the energy conservation equation. 

The different geometries have each been devel¬ 
oped to solve a specific astrophysical problem and 
hence, there are slight differences in the algorithms 
across the various coordinate systems. The first 
version of the code to be developed was the ax- 
isymmetric 2.5D version; its aim is to model the 
interaction of accreting magnetized stars with the 
surrounding accretion disk. The 2D version in po¬ 
lar coordinates was developed later on to study 
the situation where planets orbit in a magnetized 
disk. Lastly, the 3D cylindrical version of the code 
was created to study intrinsically three dimensional 
problems such as bending waves in the disk or plan¬ 
ets on inclined orbits. Each version of the code is 
based on the standard ideal MHD approach: the 
matter flow can be described with the one-fluid ap¬ 
proximation. 

In the 2D polar version of the code, the compo¬ 
nents of the velocity and magnetic field perpendic¬ 
ular to the plane of the flow are set to zero. Addi¬ 
tionally, surface density is used instead of volume 
density. The formal structure of the MHD equa¬ 
tions is the same but we make additional sugges¬ 
tions about the disk and the definitions of the (sur¬ 
face) magnetic field so that the conservation equa¬ 
tions retain the same conservative form. Our codes 
use a Riemann solver based on methods developed 


cosity terms into the code that mimic the angular 
momentum transport, but do not require high grid 
resolution. For this reason, the 2.5D version in¬ 
corporates an a-type viscosity (Shakura & Sunyaev, 
1973) which can be switched on or off depending 


by Miyoshi & Kusano (2005), modified to include 
the equation for the entropy balance. 

The 2.5D axisymmetric cylindrical version of 
the code was developed to investigate astrophys¬ 
ical processes like disk accretion, magnetospheric 
accretion onto a magnetized star and outflows 
launched from the disk and star. Such MHD flows 
may be responsible for many the observed proper¬ 
ties of young stars. In all of these processes, both 
the magnetic field of the star as well as magnetic 
fields induced in the magnetosphere of the star and 
in the accretion disk play an important role. It is 
often suggested that the angular momentum trans¬ 
port in the accretion disk is due to the magnetic 
turbulence associated with the magnetorotational 
instability (MRI) (e.g., Balbus & Hawley 1991). 
However, modelling this instability requires high 
grid resolution in the disk, which is computation¬ 
ally expensive. In many cases, it is more practical 
to allow for steady accretion by incorporating vis- 


on the phenomena being studied. 

It is less clear how the matter in the disk inter¬ 
acts with magnetic field-that is, which physical pro¬ 
cesses provide efficient magnetic diffusivity in the 
disk. Reconnection of the magnetic field lines can 
provide an efficient diffusion mechanism, but the 
process of reconnection itself also requires some 
magnetic diffusivity. Since the exact mechanism is 
poorly understood, the magnetic diffusivity is often 
parametrized in the non-ideal MHD terms. This is 
implemented in the 2.5D version of the code as a 
diffusivity module, which can also be switched on 
or off, with the diffusivity coefficient represented in 
a way analogous to the a-viscosity. 

In the 2.5D axisymmetric cylindrical version of 
the code, we split the magnetic field into a fixed 
component associated with the stellar magnetic 
field (for example a dipole field) and a field in¬ 
duced by currents in the magnetosphere and disk 
( Tanaka) , 1994). In contrast, the magnetic field in 
the 2D polar version is not split into a stellar and 
current-induced component as the code was ini¬ 
tially developed to study the interaction of a planet 
with a disk threaded by a magnetic field. In order to 
accurately compute the strength of the stellar grav¬ 
ity on the planet, the grid is allowed to co-rotate 
with the planet. This approach was suggested by 
( Kley j 1998), albeit for a different situation. 

The 3D cylindrical version of the code combines 
the two previous versions, enabling investigations 
of phenomena such as planet migration in a mag¬ 
netized accretion disk, accretion onto a magnetized 
star as well as a variety of other problems. How¬ 
ever, there is no diffusivity module implemented in 
the 3D version because the 3D instabilities which 
cause the magnetic field to diffuse are fully mod¬ 
elled. 

In this work, we describe our Godunov-type 
codes which have been implemented on several dif¬ 
ferent geometries and are designed for solving non- 
relativistic astrophysical MHD problems. Sec. [^de¬ 
scribes different Godunov approaches and the ap¬ 
proach used in our codes. In Sec. [4] we describe 
the different coordinate systems in detail. In Sec. 
[5} we describe tests of the code. Sec. [6] shows ap¬ 
plications of these codes for different astrophysical 
problems and we conclude with a summary in Sec. 
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2. The governing equations of MHD in cylindri¬ 
cal coordinates 


Here we present the full set of ideal 3D MHD 
equations. Viscosity and diffusivity are imple¬ 
mented for specific geometries and are described 


separately in Sec. 4.2.2 The continuity equation in 
conservative form is: 


dp + 1 d(rpv r ) + 1 d(pv<p) + djpv z ) _ 

dt r dr r dtp dz 


(1) 


Here r = ( r,cp,z ) are the 3D cylindrical coordinates, 
p is the density, v = (v r , v^, v z ) is the velocity, and t 
is the time. The entropy conservation equation is: 


is the net external force from the central star and 
planet. Lastly, the induction equations are: 
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This system of equations describes the evolution of 
a magnetized plasma in cylindrical coordinates. 


3. Godunov type schemes 


dips) 1 dirpsv r ) 1 dipsvcj,) <9(psv,) 

dt r dr r dtp dz 


( 2 ) 


where s = p/p r is the entropy per unit mass and y 
is the adiabatic index. The momentum equations 
in the r, tp and z directions are 
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where B = ( B r , B<p, B z ) is the magnetic field vector, 
Q = p + B 2 /8n is the total pressure, and g = -V4> 


In order to numerically solve the ideal MHD 
equations (Eqs. [lj6]), we use the Godunov method. 
In Godunov-type schemes, the central problem is to 
develop an exact or approximate algorithm to solve 
the Riemann problem. The construction of an exact 
solution to the Riemann problem is possible, but it 
can be numerically complex and is often compu¬ 
tationally expensive. The procedure of exact solu¬ 
tions of the Riemann problem has been realized by 
Ryu & Jones (1995 ) and these exact solutions are 
used for the testing of our approximate Riemann 
solver. 

Consider the one-dimensional hyperbolic system 
of equations in conservative form: 


m dTCU) 

dt dx 


(9) 


Here, U = \U\,..., U n \ is the vector of conservative 
variables and Ti'U) = {F\,...,F n } is the flux func¬ 
tion which maps each conservative variable to its 
corresponding flux. Eq. [9] can be rewritten in the 
integral form of the conservation laws 


f Uejodt = f u\ t=0 dt + r L -r R , 

J Af JAf 


( 10 ) 


where 1f ex (^,H L ,U R ) represents the exact solution 
of the Riemann problem between the left and right 
states Ul and Ur, % = x/t is a self-similar distance¬ 
time variable; is the interval along the variable 
p at which all of the waves are localized; and Tl = 
Ti'Ui ) and Tr = T ('Ur) are the fluxes for the left 
and right states, respectively. 
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Figure 1: The wave propagation diagram for the HLL method in 
the ( x, f)-plane. 


Figure 2: The wave propagation diagram for the FILLD method 
in the ( x , p-plane. 


The Riemann problem can be approximated by 
a number of discontinuities, each of which satisfies 
the Rankine-Hugoniot jump conditions: 


More advanced approximate Riemann solvers in¬ 


clude two or more intermediate states ( Gurski 
2004} [SI [2005] |Miyoshi & Kusano| [2005] ) . 


x{u 1 -u,) = r2-n. 


(id 


Here, x is velocity of the discontinuity; U\ .Ui are 
the conservative variables on either side of the dis¬ 
continuity; and T \, Ti are the fluxes of the conser¬ 
vative variables across the discontinuity. In the gen¬ 
eral case, T\ ± T(U\) and Ti + T(Un), excluding 
Ti. = T('Ul), Tr = T('Ur). In this case, the condi¬ 
tion for self-consistency of the conservation laws is 


evidently satisfied, while Eq. 12 for the calculation 
of fluxes gives T — T* , where T* is the flux in the 
interval which includes the point x — 0. If the point 
x - 0 is located at the boundary between intervals 
(that is, at the discontinuity where x, = 0), then the 
flux can be calculated for any two adjacent states. 
These fluxes are the same as the ones present in the 
Rankine-Hugoniot jump conditions. 

The approximate solution of the Riemann prob¬ 
lem is represented by the vector U(£) which ap¬ 
proximates the same conservation laws in integral 
form. The flux between the left and right states is 
approximated by 

J f-»oo 

nU£) - U L )dJ 

o 

= Tr~ 


r*oo 

am - 'U R )dt. 
Jo 


( 12 ) 

One of the simplest approximate Riemann solvers 
was proposed by Harten et al. ( 1983| ) and has 
two discontinuities separating three homogeneous 
states (as shown in Fig. [l]). One discontinuity prop¬ 
agates to the left with velocity b L , while the other 
propagates to the right with velocity b R . The ap¬ 
proximate solution is assumed to be homogeneous 
between these discontinuities: 'ZY(f) = U*. 


3.1. The Miyoshi & Kusano HLLD Solver 

Here, we describe the construction of our Go¬ 
dunov method for the equations of ideal MHD. 
As noted earlier, we use the entropy conservation 
equation instead of the energy equation. The con¬ 
servative variables and fluxes in Eq. [9] take on the 
form: 


U ={p, PS, PV X , pVy, p\’ Z , By, B z ), 

, B 2 y + B\ 
T ={pv x , psv x , pv x + p + 


pv y v x - 


B;{)By 


pv z v x - 


47T 

B x qB z 


4 n " ' 4 n 

v x B y - v y B x o, v x B z - v-B xi) \ 


Li (2005) proposed a modification of the HLLC 
approximate Riemann solver for the equations of 
MHD in which the transverse v, c-components of the 
velocity and magnetic field are assumed to be the 
same in both intermediate states. For this reason, 
in the limit B M —> 0, the HLLC solver does not con¬ 
vert to the HLLC algorithm with a purely transverse 
magnetic field; this is because arbitrary jumps of 
the y, z-components of velocity and magnetic field 
are still allowed at the contact discontinuity (as 
they are converted to the tangential component). 

The choice of entropy conservation inhibits the 
modeling of MHD flows with strong shocks at 
which entropy production occurs. Our approx¬ 
imate solution has four intermediate states—the 
wave propagation in the HLLD algorithm is shown 
in Fig. [2] The initial discontinuity between the 
states Tl / and Ur decays and the discontinuities 
appear, propagating away with velocities b L ,b R 
(the fast magnetosonic waves), a L ,a R (the Alfven 
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waves) and c (the contact discontinuity). These 
jumps separate the initial U L ,14 R and intermediate 
< U* L , < U”, ( U R , ( U R states (see Fig. [2]). The calcula¬ 
tion of the conservative variables and fluxes in the 
intermediate states is performed using the follow¬ 
ing scheme. The normal velocities across the Alfven 
and contact discontinuities are assumed to be con¬ 
tinuous 

V *xL = V 7 l. = Kr = V 7 r = C ’ 


tosonic waves: 


BmB 


V yL = V yL ■ 


xO&yL 


VxL-C 


4tt , B 2 0 

plib L -cf 


Kl = By 


g2 

pl(7> l - v xL y - —— 
4n 


pl(b L ~ cf 


Bio 

4i r 


which implies 


Analogously, we have 


( 15 ) 


Pl - PL’ Pr - Pr ’ S L - S L > S R - S R ’ 

from the Rankine-Hugoniot jump conditions. 

We write the Rankine-Hugoniot conditions for 
fast magnetosonic waves for all equations except 
the conservation equation for the jr-component of 
the momentum. For the densities in states U* L and 
'Z/jj, we obtain the relations 


'zL 
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B x qB zL 
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4n 
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4n 
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Pl ~ Pl 


b L ~ VxL 

b '- ~ Kl 


Pr = Pr 


b R - VxR 

b R - Kr 


- Pl 


- Pr 


bL ~ VxL 

b L ~c 

b R - V X R 

b R — c 


The values of v* R , v* zR , B* R , B* R are calculated using 
the same formulae, substituting the index L —> R. 
The Rankine-Hugoniot jump conditions at the 
(13) Alfven waves (with p" = Pf P* R = P * R ) can be 
solved, but only if 


From Eq. 13 and the ^-component of the momen¬ 
tum equation we obtain 


X 


Pl(v x l ~ bL ) + p R (b R - v xR ) ’ 


(14) 


where 


X = 


' 2 , B 2 B\ 

pv x + p + ^-J7 


If B 


2 \ 
*0 


PV x + P + Tn~ -47 


- b L (pv x ) L 


■ b R (pv x ) R 


Combining Eq. 13 and the Rankine-Hugoniot 
jump conditions for the entropy, we find 


S L - S L - S L > S R - S R - S R ■ 


To calculate the transverse components of the ve¬ 
locity and magnetic field in the intermediate states, 
we use the corresponding components of the mo¬ 
mentum and induction equations. In the intermedi¬ 
ate states, 1l* L and U R , these values are determined 
from the jump conditions across the fast magne- 


B x 0 B x 0 



In the first case, we adopt a L - c - \ B x p\ , w j 1 j c j 1 

A J 47 T p 1 

corresponds to an Alfven wave propagating to the 
left along the state U* L . In the second case, we take 

a R = c + -!|kL =, which corresponds to an Alfven 

f 4n pR 

wave propagating to the right along the state f U* R . 

To compute the transverse components of the ve¬ 
locity and magnetic field in the intermediate states 
r U* L *,'U R *, we use the fact that at the contact discon¬ 
tinuity 

V yl = V 7r = V v*’ V 7l = V7 = V", 

TJ** TJ** T}** TJ** TJ** T}** 

B y L ~ B yR - a y , B zL - B zK - B. , 


for B x 0 + 0. 

Plugging these into the integral conservation 
laws, we find that the transverse components of the 
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velocity and magnetic field are 



the right along U R . Additionally since a L = a R = c 
when B x o = 0, the intermediate states ( U* L *, f U R dis¬ 
appear. 

In the above derivations it is assumed that the 
inequalities 


b L < a L < c < a R < b R 


are satisfied. The case a L — c = a R , corresponding 
to /ho = 0, has been considered above. The cases 
where b L = a L or a R = b R correspond to the switch- 
on waves propagating to the left or to the right, 
assuming B tl - 0 or B, R - 0 at B x q 4 0, and are 
not considered. The wave velocities b L and b R are 
estimated as 


b L < (v x - c F ) L , b R > (v v + c F ) R , 


where cr = sign(Bvo) (Miyoshi & Kusano, 2005). 

Note that the transverse components of the mag¬ 
netic field B, are not necessarily equal {\B* L \ 4 
|Z?**| 4 \B* K \) due to the jump conditions for the 
Alfven waves in the complete system of MHD equa¬ 
tions. Here, we remove the condition T* - TifLl*) 
(in this case, for the x-component of the momen¬ 
tum equation) and hence lose the continuity of the 
magnetic pressure for the Alfven waves. 

The fluxes in the intermediate states are calcu¬ 
lated from the jump conditions 


where c f ,l, c F , R are the velocities of the fast magne- 
tosonic waves, propagating along the states 'll/, U R 
respectively. 


3.2. Correction of wave velocities 


Equations 15 and 16 show the values of the 
tangential velocities and magnetic fields in the 
intermediate states f U* L , r U* R . The numerators on 
the right-hand sides of these formulae are non¬ 
negative, while the denominators have the form 


rt = r L + b L (ui-u L ), a 7) 

rr = ri + a L cz/" - %) 

= Tl + a L CZ£* - %) + bdtr L - r U L ), (18) 
n =Tr + b R W* R - Ur), (19) 

T r = n + a R fU R - uf) 

= T r + a R W” - V* R ) + b K (U; - V R ). ( 20 ) 


One does not need to calculate all of the fluxes. 
Rather, we determine which of the states appears 
at the point x = 0 (in other words, which interval 
i-oo,b L ), ...,(b R ,oo) contains the point x = 0) and 
then calculate the corresponding flux using one of 
the formulae (jl7|)-(l20|). 

In the limit B xt) —> 0, the relationships in Eq. 15 
give the result 


V yL = VyL, 


Kl = B yL 


PL(bL - VxL) 
pl(b L - c) 2 


= B Vl 


yL 


Pl 

PL 


g2 jqI 

A = p* R (b R - c) 2 - = p R (b R - v xR )(b R - c) - . 

The formula for the left state is analogous. 

The physical sense of the last value is such that it 
should be positive because the Alfven wave should 
always propagate slower than the fast MHD discon¬ 
tinuity. The case when this denominator is zero 
corresponds to the switch-on wave, in which a fi¬ 
nite tangential field beyond the discontinuity ap¬ 
pears from the zero tangential field in front of the 
discontinuity. In the approximate Riemann solver 
it is reasonable to retain this property (the positive 
sign of the denominator) in order to avoid a non¬ 
physical change of the sign of the tangential field 
and the appearance of a singularity as A —» 0. For 
that, we require that the inequality b R > c + ^ 

yj An pR 

is satisfied, at some k > 1. In that case, we have 


Similar relationships can be derived for the z- 
component of velocity and magnetic field and for 
the fast magnetosonic wave which propagates to 


b R — c — 


k\B x0 \ 

xj^ n pR 


= b R -c- 


k\B x0 \ yjb R - c 
sj47rp R (b R - v xR ) 


> 0. 
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This can be rewritten as 


p R (b R - c)(b R - v xR ) > 


k 2 n 2 
k D 

4 n 


( 21 ) 


Comparing this relationship with that for A, we 
conclude that if the condition in Eq. 21 is satisfied 
then A > 0. Analogously: 


on the timestep by considering the wave cross¬ 
ing speeds in all cells on the grid. For our grid, 
the timesteps are typically limited by the strong 
mangetic fields in the cells near the star. For the 
calculation of the fluxes we use values from 

the intermediary timestep r l(" +l/2 , reconstructed on 
the edges of the calculated cells: 


pdb L - c)(b L - v xL ) > 


k 2 B 2 

k 

4n 


( 22 ) 


If the conditions in Eq. 21 and/or E q. |22| are not 
satisfied (after computing c using Eq. |14|) then we 
perform a correction of b R and/or b L , choosing the 
fast magnetosonic velocity to be the largest and/or 


the smallest of the roots of Eqs. 21 and/or 


Here, b R can only increase while b R can only de¬ 
crease. We then recalculate the value c, check the 


conditions in Eqs. 21 and|22[ and if necessary, once 
again apply corrections to b L and b R . 


3.3. Positivity 

We use the equation for entropy conservation in¬ 
stead of the equation of energy conservation. In 
this case, the correction of the velocities of the fast 
MHD discontinuities and the positivity of the ap¬ 
proximate solution of the Riemann solver (in the 
sense that the density and entropy are positive in 
all intermediate states) is established trivially. 


3.4. Integration in time 

To increase the accuracy of the numerical al¬ 
gorithm, we perform a two-stage Runge-Kutta in¬ 
tegration of the equations in time. In the first 
stage, we calculate the values for the intermediary 
timestep, r U" +l/2 


LI 


n+l/2 _ 


At 


=^-2A 


T£ff = FLUX CW£ 1/2 ,<W£lf)- 

The values f U L . f U R are calculated at each grid cell 
boundary using one of the previously described ap¬ 
proximate Riemann solvers. In order to increase 
the accuracy of our scheme, a slope limiter cor¬ 
rection is applied to the primitive variables at the 
edges of each cell. For the sake of brevity, we de¬ 
fine dj- 1/2 = Uj - w/ i and d,+ 1/2 = w,+i - m,-. 


u L ,i = Hi - -minmod(d/_i/2, d M p), 

hrj = Hi + ^minmod(c/ 1 _i/2,</ + ,/ 2 ). 

(23) 
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Here, the fluxes are calculated using the 

entropy-based HLLD algorithm described above. 
Formally, we can write the flux calculation proce¬ 
dure in the form - FLUX('U"/U" } ) where 

we take the left and right states in the HLLD solver 
to be U L = 'Ui and U R = U i+ \. 

In the second stage, the values for the full 
timestep, t” +1 = t" + At, are calculated as 


Here, the top indicies are dropped and we take e = 
0.5. In our algorithm, we reconstruct the primitive 
variables p, s, v x , y v , v z , B x , B v , B z . 

4. Numerical algorithms in cases of different ge¬ 
ometries 

4.1. Cartesian geometry 

In a 2D Cartesian (x, y) geometry, the MHD equa¬ 
tions take the form 


The timestep At is set by the Courant-Friedrichs- 
Levy (CFL) condition, which places an upper limit 


ffU dT dQ 
dt dx By 
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Unlike the ID case, the vector of the conservative 
variables 14 has a total of seven components. The 
sixth and seventh components are B x and B y , which 
must satisfy the divergence free condition 


dB x 8B y 
dx dy 


(24) 


The numerical algorithm for 2D is analogous to the 
ID case. The integration in time is performed using 
the two-step Runge-Kutta method 


14 


71+1/2 _ 

-'J 


At 


- ni n _C F" 

“ 2Ax i+l ' 2 ’j 

At 


-T'UnJ 




14 n+1 =14" - -(T‘ 

'■J ’J A-F i 


At 


71+1/2 


-r; 


71+1/2 


J Ax^ i+l/2j i-l/2j 


) 


At 


/ >o//+l/2 y^)Il+ 1/Z \ 


->n+l/2 


At 


the 


F LUX( 14 "j, 14 " +lj ), 


cjrn 

'*+1/2 J 


same 


time, 

Sh + 1/2 = FLUX( 14 l r 14 l j+ 1). 

An analogous algorithm for the calculation of 
fluxes is also used during the second stage of the 
time-advancement computation. 

To ensure that the magnetic field is divergence- 
free in 2D, we use the method proposed by Bal-| 
sara & Spicer (1999). In application to problems in 
Cartesian coordinates, this approach is the follow¬ 


ing: if Eq. 24 is satisfied, then one can represent 
the x and y-components of the magnetic field in 
the form 

dA z dA z 

Bx — jQ By — n 

oy ox 


At the same time, —- = -E z . When the equations 
at 

are discretized, the conservative variables 14 are 
determined in the cells while the x and y fluxes, 
T and Q, are computed on the sides of the cells. 
The corresponding components of the flux function 
vectors T and Q are 


F 6 = 0, F 7 = -E z , G 6 = E z , G-i = 0, 


where E z = -(v x B y - v y B x ) is the z component of 
the electric field (multiplied by the speed of light). 
Thus, the z component of the electric field is calcu¬ 
lated twice using two different methods: the first as 
-f 7./+I/2,/, and the second as G^ij+i/z- If 4 - is deter¬ 
mined on the nodes of the grid (i + 1/2, j + 1 /2), the 
components of the magnetic field can be calculated 
on the sides of the cells as 


B 


x, 1+1/2,;' — 


U.i+lj 


' A z ,iJ 


Ay 


B 


’y,i,j+ 1/2 - 


A zx j+1 A 


Z,l.J 


Ax 


At the same time, the divergence-free condition will 
be satisfied if we require 


£.'-■/+1/2J - B xJ -U2J + ByjJ+l/2 - ByjJ-l/2 _ 

Ax Ay 

The magnetic field in the cells are then computed 
from the calculated values on the cell boundaries 

B x ,i J =y(R.t,i-l/2J + B x ,i+l/2,j), 

By,i,j ~ ^(ByJJ- I /2 + Byij + 1/2) . (25) 


To calculate A z j+i/ 2 j+i /2 at the n + 1 th timestep, we 
determine the ’-component of the electric field at 
each node. The simple variant proposed by Balsara 


& Spicer (|1999|) consists of averaging the values 


-F-j and G(, over the cell sides surrounding the node 
at (i + 1/2, j + 1 /2). For a homogeneous grid, this is 


F BS 

c z,i+l/2j+l/2 _ 

^ {E z ,i.j +\/2 + F z i + ij + \/2 + E z i + if2J + B z j +1 /2J+ 1) — 

^(f?6,!j+l/2 + G6 ,i+1j+1/2 - F lti+ i/2J - f7,i+l/2J+l) • 

(26) 


In Gardiner & Stone fl2005p it was noted that 
this approach is inconsistent with the numerical 
integration algorithm in the case of ID problems 
(where the solution does not depend either on x or 
on y), and that the proposed procedure does not 
guarantee consistency. Gardiner & Stone (2005) 


propose a modified form (for homogenous grids) 


= 2 E 


.BS 


7 gs 

■^,1+1/2,.7+1/2 ~ ^z,i+l/2,/+l/2 

~t E Z j+ij + E z j+ij+i 


j( e m 


-z,i,j +1 


)■ (27) 


To calculate E-, we use either Eqs. 26 or|27| Eq. 27 


has applications in damping out small-scale oscilla¬ 
tions of the magnetic field. 

Thus, at each stage of numerical integration of 
the MHD equations, the magnetic field in the (x,y) 
plane is calculated along the following algorithm: 

1. Use an approximate Riemann solver to find the 
z-components of the electric field at the sides 
of cells -F-JJ+ 1 / 2 J, G( iX j + 1 / 2 * 

2. Use Eq. 26 or 27 to calculate the electric field 
in the nodes of the grid £ : ,!+i/2j+i/2- 

3. Calculate the z-component of the magnetic 
field potential in grid nodes, A z j + i/ 2 j+i/ 2 - 

4. Calculate B x , i+ i/ 2 j and B yXj+ 1 /2 . 

5. Using Eq. 25 calculate B xx j and B y jj. 
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4.2. Axisymmetric 2.5D cylindrical geometry 


To solve the 2.5D MHD equations (with mag¬ 
netic diffusivity and viscosity) in the disk, we split 
the different physical processes into several differ¬ 
ent modules. These are: (1) the “ideal MHD” 
module in which we calculate the dynamics of 
the plasma and magnetic field with dissipative 
processes switched off; the “diffusivity” modules 
(2) and (3) which calculate the diffusion of the 
poloidal and azimuthal components of the mag¬ 
netic field for frozen values of the plasma velocity 
and thermodynamic parameters (density and pres¬ 
sure); and the “viscosity” module (4) which com¬ 
putes the viscous dissipation in the disk. 

(1) In the hydrodynamic module, the ideal 
MHD equations are integrated numerically using 
an explicit conservative Godunov-type numerical 
scheme. In our numerical code, the dynamical vari¬ 
ables are determined in the cells, while the vector- 
potential of the magnetic field, A#, is determined 
on the nodes. To calculate the fluxes between the 
cells, we use the previously described approximate 
Riemann solver. 

For better spatial resolution, we perform the 
reconstruction of the primitive variables to the 
boundaries between the calculated cells. We use 


two types of limiters: a) the MINMOD limiter (Roe 


1986]), and b) a limiter which reconstructs the func 
tion with 3rd order accuracy (which is similar to 
the limiter developed by Koren 1993). Integration 
of the equations in time is performed with a two- 
step Runge-Kutta method. 

To guarantee a divergence-free field in the 2.5D 
version of the code, the ^-component of the mag¬ 
netic field vector-potential A^ (or magnetic flux 
function *F = rA,,,) is calculated at each timestep 
and then used to obtain the poloidal components 
of the magnetic field ( B r , B z ) in a divergence-free 
form ( Toth| [2000p . In other words, the condition 
V • B = 0 is guaranteed to be satisfied to within ma¬ 
chine accuracy. The calculation of A,,, is performed 
using either the algorithm proposed by Bals ara &| 
Spicer ( 1999|) or the modified form proposed by 
Gardiner & Stone (2005). 

(2) In the module where the diffusion of the 
poloidal magnetic field is calculated, we numeri¬ 
cally integrate the equation for the ^-component 
of magnetic field vector-potential (or magnetic flux 
function V F), using an explicit approximation. 

(3) In the module where the diffusion of the az¬ 
imuthal component of the magnetic field and (4) 


the viscous stresses are calculated, we add diffu¬ 
sive terms which approximate the corresponding 
spatial derivatives to the fluxes calculated by the 
ideal MHD module. 

The ideal and diffusivity modules are verified 
separately using a number of standard tests. The 
setup and results of these tests are discussed in Sec. 

m 


4.2.1. Ideal MHD module 

The ideal MHD equations in the axisymmetric 
cylindrical geometry are derived from the full set 
of 3D MHD equations by requring that d/drp - 0 
(see Sec. 4.31. The equations are written in a form 
which does not have singularities on the right hand 
side. For the r- and ^-components of the equation 
of motion and ^-component of the induction equa¬ 


tion, we substitute the operator r, representing 

r or 

divergence operator in cylindrical coordinates. As 
a result, the ^-component of the equations of mo¬ 
tion takes on the form of angular momentum con¬ 
servation while the ^-component of the induction 
equation takes on its typical form. 

In discretizing the cylindrical MHD equations, we 
take into account the fact that near the symmetry 
axis (where r —» 0), different variables have differ¬ 
ent asymptotic behavior. For example, the values 
p,p,v z ,B z —» const as r —> 0 whereas the values 
v r ,v$,B r , Bq —> 0 scale with r. Analogous behav¬ 
ior is shown by the fluxes corresponding to these 
variables. In particular, in the approximation of the 
^-component of the equation of motion (the angu¬ 
lar momentum equation), we multiply the fourth 
equation (Eq. [ 4 ]) by r 2 3 and integrate it along r and 
z in the limits of the calculated cell (i, j ). The equa¬ 
tion is then rewritten in the form = a LJ r and 


B B 

G = pv^Vr - 1 ^- = Ayr, where ay and Ay weakly 

depend on r. Here we have dropped the indices on 
F and G. 

We have 


d 

dt 




ijF dr + A r (-) + A z 




rdr = 0. 


Next 


f 


dijrdr = 


2 ,2 

/ \ 'i-1 /2 + *7+1/2 f i'-l/2 + 7*1+1/2 , v 

\ a ui) -x-x- (n + 1/2 - n-i/2) = 


v 4>Ai (r 2 ),. A n , 
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where 


: ( a i-j ) 


^ 1 — 1/2 + r,+i /2 


2 , 2 
r i—1/2 + r i+l/2 


n- 


Ar,- =r i+ i / 2 - r,-_i/ 2 . 


(28) 


Analogously, 


J' A Lj± i, 2 r 3 dr = G iJ±xj2 (r 2 ). Ar,-. 

Thus, the angular momentum conservation equa¬ 
tion for the cell (i, j) has the form: 


d(pv<f>)i,j rf +U2 F i+ i/ 2 j - r^-ii 2 F i-il 2 ,j 


dt 


(r 2 ) ; Ar,- 

G/j+1/2 - G/j-1/2 


AZi 


= 0. 


(29) 


It is clear that this approximation of the operator 
1 d 

— — r 1 and more precisely of the multiplier (r 2 ), 
r 2 dr ' H 

significantly differs from analogous terms (for ex¬ 
ample, of r 2 ) only in the vicinity of the symmetry 
axis. 

The flux densities Fi±t/ 2 ,j and G iJ± 1 / 2 , which are 
used in Eq. [4j are calculated using the ID algo¬ 
rithm for the approximate solution of the Riemann 
problem. 


4.2.2. Viscosity and diffusivity 

Here we discuss the treatment of viscosity and 
diffusivity. We use equations of motion in cylindri¬ 
cal coordinates in the form: 


dpv r dT rr dT rz 


Tthtb ~ T„ 


dt 


dr 


dz 


+ Pgr, 


dpvtf, 1 dr~T r( p dT zl p 

+ — —-—- + —— - 0 , 


dt r 2 dr 


dz 


(30) 


dpv z 1 drT rz cYT, z 

+ “ —— + = Pgz- 


dt r dr dz 
The stress tensor = T + r,/ f consists of an ideal 
part T ik : 

( B 2 \ B,B k 

Tik = pViV k + I p + ^ I o ik - ~ , i,k = r,(f>,z , 

which is used to derive the momentum equations 
in Eqs. [3] [4] and [5] and a viscous part r,* which 


takes into account small-scale turbulent velocity 
and magnetic field fluctuations. 

We assume that the stress due to the turbulent 
fluctuations can be represented in the same way as 
the collisional viscosity by substituting in the tur¬ 
bulent viscosity coefficient. The components of the 
tensor r* in axisymmetric cylindrical coordinates 
are: 


da> dot v r 

T r<j> = ~pr— ,r , 0 = -pr— ,t h = - 2 p— , 
dr dz r 

dv r dv r 

Trz = ~Pyr ’ T 't = ’ 

dz dr 


(31) 


Here to = v (!) /r is the angular velocity of the plasma 
and p is the dynamical turbulent viscosity. In 3D, 
we only consider the rtf> and zip components of the 
stress tensor as these are expected to be dominant; 
similarly, in 2D polar coordinates, we only consider 
the rip component. 

Separating out the viscous stress in the <p compo¬ 
nent of the momentum equation gives 


dipvcj,) d(r 3 T r< fi) 


1 <9(sin 2 GT Z $) 


dt 


1 d 


r 3 dr 


da) 


AZ W sin0 HT + 


dr 


r sin 2 6 
\ _ d_ 

in 2 e se 


d6 


da) 


v pS ,n#-j. 

(32) 


where T r(j) and T 7ll< are components of the inviscid 
part of the stress tensor. 

We assume that the plasma has a finite magnetic 
diffusivity. That is, the matter may diffuse across 
the field lines. We also assume that the finite diffu¬ 
sivity of the plasma is due to the small-scale turbu¬ 
lent fluctuations of the velocity and the magnetic 
field. The induction equation averaged over the 
small-scale fluctuations has the form 

hD 

-V x (v x B) + V x E 1 ' = 0 . (33) 

dt 


Here, v and B are the averaged velocity and mag¬ 
netic fields, and E f = - (v' xB') is electromotive 
force connected with the fluctuating fields. 

Since the turbulent electromotive force E 1 is con¬ 
nected with small-scale fluctuations, it is reason¬ 
able to suppose that it has a simple relation to 
the ordered magnetic field B. If we neglect the 
magnetic dynamo cr- effect ( Moffat} 1978), then 
(v' x B') = - 77 V x B, where r\ is the coefficient of 
turbulent magnetic diffusivity. Equation (331 now 
takes the form 
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- Vx(vxB) + Vx ( 77 V xB) = 0. 

dt 


(34) 



























We should note that the term for FJ formally coin¬ 
cides with Ohm’s law 

c c 2 f 

J = —VxB = -E f . 

An Anri 

The coefficient of turbulent electric conductivity 
cr = c 2 / Am). 

To calculate the evolution of the poloidal mag¬ 
netic field it is useful to calculate the ^-component 
of the vector-potential A or magnetic flux function 
'F = rA<p. Due to the assumed axisymmetry, 


1 W _ 1 d'F 
r dz’ ' r dr 


(35) 


Substituting B = V x A into the induction equa¬ 
tion gives the equation for the 0 component of the 
vector-potential 


d'V Id Iff 3 2 v F\ 

dt ^ \ dr r dr dz 2 / 


= r(\ x B ) 0 . 


(36) 


The azimuthal component of the induction equa¬ 
tion gives 


dlii, d IrjdrBA d I dBA 
dt dr\r dr ) + dzV dz ) 


(V x (v x B )) 0 . 

(37) 


The viscosity and magnetic diffusivity lead to dis¬ 
sipation of the kinetic and magnetic energies, con¬ 
version into thermal energy and a corresponding 
increase of the entropy. In our simulations we have 
neglected viscous and Ohmic heating. We have also 
neglected radiative cooling. Thus, the main effect 
of the viscous terms is the transport of angular mo¬ 
mentum outward which allows matter to accrete 
inward to the disk-magnetosphere boundary. The 
main effect of the diffusion terms is the transport 
of matter through magnetic field lines. 

The dissipative terms t Z(! , and, IT are incorpo¬ 
rated into the numerical algorithm by way of a flux 
correction in the finite-difference equations. The 
modified fluxes F 4 , G 4 are included in the angular 
momentum conservation equations, and the modi¬ 
fied fluxes Ft, Gt (components of the poloidal elec¬ 
tric field) are included into the ^-component of the 
induction equation. The approximation of the vis¬ 
cous stresses takes the form: 

p _ pid ,, .. a>i+1 J ~ 

^ 4,i+l/2j - *4,i+l/2,j ~ Bi+t/2,j r i+l/2 -T--. 

Ar i+ 1/2 

^ d .. OJij+i-OJij 

G 4 ,ij+l/2 - G 4,i,j+l/2 -BiJ+l/in ^ ■ 


The approximation of the poloidal electric field has 
the form: 


F l,i+l/2,j ~ 
G 7,U+l/2 = 


^ 7,i+1/2 ,j 1i+W 


G l,i,j+l/2 tli,j+l/2 r i 


t"i +1 B t 1 j V\B 


(rAr),-+i / 2 

d(pjj +1 — B 
A+/+ 1/2 


AJ 


A ,j 


Here, F 14+1/2 j and G 1JJ+ i/2 are the z- 
and r-componenls of electric field, and 
F‘ii+ 1 /2 p G ‘ii j+ 1 p are ^e same components cal¬ 
culated in the approximation of the ideal plasma 
conductivity using the algorithm described above, 
(rAr) i+ i /2 = (r 2 +1 - r 2 )/ 2 . 

The coefficients of the dynamical viscosity at the 
boundaries between calculated cells are calculated 

according to Bi+i/2,j = . Here, and 

Bi,j "t Bi+lJ 

Hi+\j are the coefficients of the dynamical viscosity 
in the grid cells. Analogously, we calculate Hij+ 1/2 
and the coefficient of the dynamical viscosity at the 
boundaries between cells. 

The azimuthal component of the electric field, 
responsible for the evolution of the magnetic 
flux function 'F, is approximated using a special 
method. The azimuthal component of the turbu¬ 
lent electric field E ^ is presented in the form 


I d 1 dW d 2 ^) 
T d-rrd^+d^j 


and approximated by 


E* 


0,/+1/2,7+1/2 ~ 7 7 i+ 1/2J+1/2 


n+ 1/2 


Ar ;+ i /2 

'f < (+3/2j+l/2 - 'f , i+l/2j+l/2 

(rAr) M 

^ 1 + 1 / 2 , 7 + 1/2 “ 'Pi- 1 / 2 , 7 + 1/2 
(rAr), 

1 f'Pi+l/2,7+3/2 -'f'i+l/2,7+I/2 


A+ 7 + 1/2 \ 


A+7+1 


V F, 


i+1/2.7+1/2 ■ 


V F, 


i+1/2,7-1/2 


A Zj 


(38) 


A corresponding term is added to Eq. 36 which 
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takes on the form 


iim +1 

*+l/2j+l/2 


VT/n 

T /+ 1 / 2 j+ 1/2 


- ?7/+1/2,y+1/2 


n+ 1/2 


Ar LAr/+i/2 

'f'i+3/2j+l/2 - 'P;+l/2j+l/2 

(rAr) i+ i 
*y;+i/2j+i/2 - *y;~i/2j+i/2 \ 
(rAr), } 

1 /'Pi+l/2j+3/2 - 'P;+l/2j+l/2 

Azj+i/ 2 \ Az J+ i 

'f , i+l/2j+l/2 - 'f , i+l/2J-l/2^ 


After calculating the electric field components, 
the magnetic field components at the sides of cells 
(for example, B r at (i + 1/2 ,j,k)) are calculated us¬ 
ing the formula 

D/7+1 _ D/7 -pn _ EV7 

n r,i+l/2,j,k n r,i+l/2,j,k ^ z,i+\l2,j+\/2,k z,/+l/2J-l/2,£ 


At r i+ i /2 A(f) 

pn _ pn 

^cp, i +1 /2,j,k- 1/2 c 0,!+1 12 ,j,k +1 /2 

Ac 


= 0 . 


At; 


-h+i/ 2 ^. 


0,/+1 /2,j -\-1 /2" 

(39) 


The values of E'£ i+l/2 . +1/2 on the right hand side are 


calculated using Eqs. 26 and 27 At the first step 
of the time integration, we use a time-step of Af/2 
and calculate the values of't'"/. 1 /, 2 . , n . 

1 + 1 / 2 , 7 + 1/2 

4.3. 3D cylindrical geometry 

We have also developed a fully three dimensional 
version of the code in cylindrical coordinates. To 
calculate the fluxes between cells, we use the HLLD 


The other components B^-j+i j 2 ,k and B z jj^+ 1/2 are 
calculated in a similar fashion. This ensures that 
the magnetic field (B r j +U2J , k , %, ; -+i/ 2 ,*, B zXjM1/2 ) is 
divergence-free (so long as the field is initially 
divergence-free at the beginning of the simulation). 
As a next step, the magnetic field in the cells is cal¬ 
culated from the values on the cell boundaries: 

B rAjJc = \(Br Jr . 1 l2J* + B rj+ll2J t), 

B <P jjy=-(B M - l/2 ,k + B^ j+l/2ik ), 

B z ,i,j,k — y ( BzJ. j.k -1 /2 + B-.i.j.k +1 /2^ ■ 


method as in the 2.5D code (Miyoshi & Kusano 


2005). To correct the fluxes, we use the MINMOD ^.4. 2D polar geometry 


flux corrector, which is applied on the primitive 
variables p, s, v r , vy, v-, B,., B t ,„ B,. To integrate the 
equations in time, we use the two-step Runge-Kutta 
method. 

The main differences between the 2D and 3D 
codes lie in the algorithm which ensures that the 
magnetic field is divergence-free. Namely, the 
divergence-free condition is guaranteed by calcu¬ 
lating the components of the magnetic field on the 
sides of cells from the induction equation: 

Y + V x E + 0. 
at 

The components of the electric field are calculated 
at the edges of the cell as in the axisymmetric case. 
For example: 


To describe the interaction of the accretion disk 
with a planet in the presence of a magnetic field, 
we use a model which describes the processes 
in terms of surface variables (density and pres¬ 
sure). The equations for this model are obtained 
by integrating the three-dimensional MHD equa¬ 
tions along the ’-direction which is perpendicular 
to the disk plane. For simplicity, we assume that 
the --components of the velocity and the magnetic 
field are zero and there is no field outside the disk. 

The governing equations are derived from the 
ideal MHD equations (Eq. [ 6 ]) by integrating in the 
vertical direction: 


d? 1 dr?, iv 1 d?v 0 
dt r dr r dcf> 


= 0 , 


(40) 


-Z,i+ll2,j+l/2,k 


1 / 

T [E z ,iJ+l/2,k + E z j+lj+l/2,k 


+ Ez,i+l/2,j,k + E z j +l /2j+l,k) , 

d? iv d 9 

dt + dr Sl V +n + 

where the values on the right hand side are the 
corresponding fluxes calculated by the Riemann 
solver: £/,;,;+ 1 / 2 , i is determined between cells (i,j,k) 
and i,j+ 1, k and so on. 

V 

1 d I b r bA 

+rj = 


2 f 


8 tt 


b 

An 




bl-bl 

<P r 

4 n 


-? 


GM t 
r 2 

(41) 


+ Z/r, 
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dt 

\_d_ 

r d(f> 


IM 

r 2 dr \ 




Svt; + n + 


v 


2v,.v 0 - 

8 n 


b r bA 
4 :t ) 


= *U- 


dIS 1 drI,Sv r 1 52Sv 0 

dt r dr r d(f> 


(42) 

(43) 


Here 2 = f pdz is the surface density; II = f pdz 
is the surface pressure; b = (b r ,b<p,b z ) is the z- 
averaged analog of the mangetic field; S = n/ 2 r is 
an analog of the entropy; 2 / r , and 2 / (t are the com¬ 
ponents of the force per unit surface area acting 
from the planet on the disk matter. These resemble 
the standard form of the 2D equations where the 
z-derivative in the 3D equations is set to 0. 

The induction equations take on the form 


dy ^' b ‘ +~ r ^ VM!■>£,. - v r b^) = 0, (44) 

^ \/Mv,.b 0 - vpbr ) = 0. (45) 

Here, h„, is the “magnetic thickness” of the disk, 
which is the characteristic scale for the magnetic 
field distribution in z-direction. This scale is not 
identical to the disk height h which characterizes 
the z-distribution of the density and pressure in 
the disk. We take h m to be the same across both 
equations. Since the magnetic field in the disk is 
unlikely to change sign in the vertical direction, the 
approximation seems to be reasonable. The next 
suggestion which helps to close the system of equa¬ 
tions is in the dependence of the magnetic thick¬ 
ness of the disk h m (r) on the radius. In particu¬ 
lar, if we assume h m is constant, then the induction 
equations will also take the standard form (in the 
same sense that J1 =0). To calculate the magnetic 
field and volumetric density and pressure from the 
z-integrated values, we assume a disk thickness at 
some characteristic radius and back out the volu¬ 
metric terms. 


5. Tests of the Code 

We have performed multiple tests of the code 
across different grid geometries. Below we show 
two examples of such tests. In the first, we test the 
ideal MHD module (with the viscosity and diffusiv- 
ity modules switched off). In the second, we test 
the diffusion module separately. 



Figure 3: The test of the ideal MHD module of the code with the 
“rotor problem” at different grid resolutions. Top Panels: Density 
distribution (color background) and field lines (solid lines) with 
grid resolution 100 x 100 (left panel), 200 x 200 (middle panel) 
and 400 x 400 (right panel). Bottom Panels: The density contour 
lines. 


5.1. Test of the ideal MHD module 

5.1.1. The 2D rotor problem 

The first test is the standard “rotor problem” in 
Cartesian coordinates. This test has been used by 
a number of authors for testing MHD solvers which 
use the energy equation (Balsara & Spicer, 1999; 
Ttathl [2000] ). We use this test to check the ideal 


MHD module of our code (with viscosity and diffu- 
sivity switched off). 

The ideal MHD equations are solved on a reg¬ 
ular Cartesian grid in the region -0.5 < x < 0.5, 
-0.5 < y < 0.5 with grid sizes Ax = Ay = 1 /N, where 
N = 100,200,400 for the three different tests. At the 
beginning of the simulations, t = 0, the pressure in 
the region is constant, p = 1 , and the magnetic field 
is homogeneous, B x = 0, If -5. In the center, there 
is a circle of higher-density matter (po = 10) with 
radius r 0 = 0.1 (where the radius is r = a/x 2 + y 2 ). 
The matter in the inner circle initially rotates as 
a solid body with angular velocity too = 20. For 
r > n = 0.115, the density is p\ = 1 and the matter 
is at rest. In the ring in between these two regions 
r 0 < r < r\, the density and velocity are linearly in¬ 
terpolated between the values at r = r 0 and r - r\. 

The equations of ideal adiabatic MHD are 
solved with the previously described Godunov-type 
scheme. The timesteps are calculated from the 
condition At = 0.4A i C fl where At chL is the maxi¬ 
mum timestep allowed by the Courant-Friedrichs- 
Levy (CFL) condition. The results of the simula¬ 
tions at t = 0.15 are shown in Fig. [3] Evidently, the 
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Figure 4: Slices of the density p, pressure p, kinetic energy, Ek 
and magnetic energy E,„ along a horizontal line at the center of 
the rotor at time t = 0.15. 


Figure 5: The test of the ideal MHD module with the Low-type 
solution using a homogeneous grid with N r , N z = 80 x 160. From 
left to right: contours of the pressure, density and magnetic flux 
functions at the moment of time when the scaling parameter is 
a = 6. In all of the plots, the solid line shows the numerical 
solution while the dashed line shows the analytical Low-type 
solution. 


P p ¥ 



Figure 6: The same as in Fig. [5] but using an inhomogeneous 
grid with N r , N z = 40 x BO. 


density and the field line distribution is very similar 
in all three of our cases. The bottom panels of Fig. 
[3] show selected streamlines with numbers which 
clearly show the similarity of the results across the 
various grid sizes. Fig. [ 4 ] shows slices along the 
y-axis for the density, pressure, kinetic energy, and 
magnetic energy for each of the three grid sizes in 
addition to a high resolution grid with N = 800. At 
the highest resolutions, the models are nearly in¬ 
distinguishable and the convergence of the results 
is evident. 


5.1.2. 2D Low-type analytical solutions 

To test the ideal MHD module, we also compare 
our numerical calculations with analytic solutions 
found by Low (1984). One of them has been used 


in 


Stone & Norman (1992) for testing the corre¬ 


sponding module of the ZEUS code. It is conve¬ 
nient to write these solutions in spherical coordi¬ 
nates (as we will do below), and then to convert 
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The magnetic flux function has the form 



V = Cf{f - ^) 2 sin 2 6 . 

The poloidal components of the magnetic field are 
calculated along the known formulae 


Figure 7: The test of the ideal MHD module using the Low-type 
solution on a homogeneous grid with N r ,N z = 80 x 160. The de¬ 
pendence of the radial component of velocities and components 
perpendicular to the radius on the spherical radius for all grid 
cells. 



Figure 8: The test of the ideal MHD module with solution of 
|Stone & Norman ( [1992} on homogeneous grid with N r ,N : = 
70 x 140. From left to right: pressure, density and magnetic flux 
function at the moment of time when the scaling parameter is 
a = 6. In all plots, the solid line shows the numerical solution, 
dashed line - analytical solution. 


them to cylindrical coordinates. 

These solutions are self-similar and assume that 
the magnetic flux function has a dependence on 
the time and radial coordinate of the form £ = 
R/a(t). Here a is a time-dependent scaling param¬ 
eter and £' is the self-similarity radius. The ve¬ 
locity has only a radial component (in the spheri¬ 
cal coordinate system) which depends linearly on 
the radius. The function a(t) satisfies the equation 
1 / dci \ ^ ex 

- — h— - P, and the magnetic flux function 


2 \ dt 

satisfies the Grad-Shafranov equation. The other 
variables are connected with a(t ) and 'F and their 
derivatives by algebraic relationships. The first so¬ 
lution proposed by Lou (1998) and used for testing 
does not have an azimuthal component of the mag¬ 
netic field. The second solution used in iStone &i 


Norman (1992) has a non-zero azimuthal magnetic 


field but the motion of the plasma is inertial. 

In the first solution, Lou (1998) chose growing 
solutions of the function a(t) 


da 

dt 



(46) 


2Ccos0, , 

b r = — -(f-tir. 


B ff = - 


a 

2 C sin 6 




The initial pressure and density are 


P = 


C 2 


(5? 0 -^ 2 )(?-? 0 )? S m 2 e-e 0 

2 C 2 % 


na 2 (a^ + GM/p 2 ) 


7(r-£o)Tsin z 0 + ^ ■ 


And the radial velocity is given by: v R = 

*#-!)• 

We take the following parameters for the prob¬ 
lem: GM = 1, C = 1, or = 1, p = 1, £ 0 = 6, the 
adiabatic index is y — 4/3. 

The solution is numerically computed over the 
region 0 < r < 8, -8 < z < 8 with a rectangu¬ 
lar region, 0 < r < 1,-1 <z< 1 excised from 
the simulation. We use two grids: a homogeneous 
grid, with N r x N z = 80 x 160, and an inhomoge¬ 
neous grid, 40 x 80. In both cases, the excised rect¬ 
angular region is homogeneous, with dimensions 
10 x 10, though values are not computed in this re¬ 
gion. Thus, the cell sizes of the homogenous grid 
are Ar = 0.1 and Az = 0.1. For the inhomogeneous 
grid, the grid size is increased in a geometrical pro¬ 
gression with the power q = 1.05. The integration 
time-step has been chosen automatically such that 
At = 0.5A t C FL- The function a(t) is determined by 
numerically integrating Eq. [46] with the initial con¬ 
dition a( 0) = 2 at sufficiently fine grid resolution. As 
an initial condition, the analytical solution is used 
with the parameter a = 2. For boundary conditions, 
we also use the analytical solution for all variables 
with the parameter a corresponding to the moment 
of time. The simulations were performed up to the 
moment of time when a = 6. Fig. [5]and Fig. [6]from 
left to right show the pressure, density and mag¬ 
netic flux function at the moment of time when the 
scaling parameter is a = 6. Fig. [5] shows the result 
of the simulations on a homogeneous grid, while 
Fig. [6] shows the results for the inhomogeneous 
grid. In all of the plots, the solid line shows the 
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The density and pressure are given by 



Figure 9: The test of the ideal MHD module using the solution 
of |Stone & N orman (19921 on homogeneous grid with N r ,N : = 
70 x 140. The dependence of the radial component of velocities 
and components perpendicular to the radius on the spherical 
radius for all grid cells. 


numerical solution and the dashed line shows the 
analytical solution. Fig. [7] shows the radial and 
transverse components of the velocities as a func¬ 
tion of the spherical radius. One can see that the 
velocity is directed along the radial direction and is 
proportional to the radius with the high accuracy. 

In the second solution flStone & Norman] 1992 ) 
the scaling factor linearly depends on time, a(t) = 
t ^frj. In this case, the magnetic flux function has 
the form 


*F = 


A 0 | po+ 
0 


sind^ 

~w 


cos /Id sin" 6 


t<Zc 


The poloidal components of the magnetic field (in 
spherical coordinates) can be calculated using the 
analytic solution 


Br = 


2Aq cos 0 I 

~1^[ P0 + 


sindd 

~w 


- COS /Id 


0 




Ba = • 


/IAo sin 6 
ciR 


sin dd + 


cos dd sin dd 
dd ” d 2 d 2 


0 


f >fc 


The azimuthal component of the magnetic field is 

d'P 


dAo sin 6 
aR 


B(h — 

a z R sin 6 

sin dd \ 

Po + ’ - cos dd d < %c 




(47) 




p = p s + 


A 2 0 Po 
2nGMR 3 


(4-d 2 d 2 )|p 0 + 


sindd 


- cos Ag sin 6 , 


Ps — < 


/GM\ 3 




Of.- 


A oPO, 


P = Fv + 73tj( 2 - d 2 d 2 ) I Po + r: ^ L - COS dd | sin" 0, 


4tiR 4 


sindd 





GMp s 



AR 


Ps = ■ 

IdopRl 

i2GM IR {) ^\ 


6 u 4 ^ 

{ 3^o 1 f / / 


do = 10uHp exp 


2GM) 

3^n i 


, w p = 1.673 x 10“ 2 V 


The radial velocity is simply \\ — d = /?/f. 

The parameters of the problem are assigned the 
following values: M — 2 x 10 33 g, 77 = 5.24 x 
10 _8 sec" 2 , d = 5.54 x 10 -11 cm -1 , dr = 1-104 x 
10 n cm, v = 2.42 x 10 20 cm 3 sec _2 g^ 1/3 , Ao = 1.5 x 

10 21 Gcm 2 , p 0 = cosddc - Sm = 1.01327, R 0 = 

ddc 

10 n cm, the adiabatic index is y = 4/3. 

The equations of ideal MHD are integrated in the 
region 0 < r < 3.5 x 10 n cm, -3.5 x 10 n cm < z < 
3.5x 10 n cm, where an inner rectangle with size 0 < 
r < 10 H cm, -10 n cm < z < 10 n cm has been excised 
from the region. We use the homogeneous grid 70 x 
140, and the grid in the excised region is 20 x 20. 
Thus, the cell size is A r = 0.05 x 10 n cm and Az = 
0.05xl0 n cm. The time-step is chosen automatically 
such that At = 0.5Af C rx- As an initial condition, we 
take the analytical solution described above for a = 
2. We use the analytical solution for the boundary 
conditions, taking the parameter a to correspond to 
the moment of time. The calculations are done up 
to the moment when the scaling parameter reaches 
the value a = 6. As before, Fig. [8] and Fig. [9] show 
the results of simulations at the moment when the 
scaling parameter becomes equal to a = 6. 


5.1.3. The 3D rotor problem 

As with the 2D code, we test the 3D code with 
the rotor problem. The 3D MHD equations are 
solved numerically in the cylindrical region r < 0.5, 
-0.5 < z < 0.5 and the initial conditions are deter¬ 
mined using the same initial conditions as in the 
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Figure 10: Test of the ideal MHD module of the 3D version of 
code with the “rotor problem” on a grid with N r ,N^,N z = 500 x 
500x20. Left Panel: Density distribution (color background) and 
field lines (solid lines). Right Panel: The pressure distribution 
(color background) and field lines (solid lines). 


2D test. In the beginning of the simulations, t = 0, 
the pressure in the region is constant, p = 1 , and 
the magnetic field is homogeneous, B x = 0, B y = 5. 
In the center there is a circle with radius ro =0.1 
where the matter density is po = 10 and the mat¬ 
ter rotates as a solid body with angular velocity 
w 0 = 20. At r > r\ = 0.115, the density is p\ = 1 and 
the matter is at rest. In the ring r 0 < r < r\, the den¬ 
sity and velocity are linearly interpolated between 
those at r = ro and r = r\. 


In the test we use a homogeneous grid with 
N r ,N l j l ,N z = 500 x 500 x 20. The integration time- 
step is chosen automatically in such way that At - 
OAAtcFL■ The results of the simulations at time 
t = 0.12 are shown in Fig. 10 One can see that 
the density, pressure and the field line distributions 
are very similar to 2D case. 


5.1.4. 3D Low-type analytical solution 

To test the ideal MHD module, we also use the 


analytical solution described in Sec. 5.1.2 Cal¬ 
culations were performed in the cylindrical region 
r <6. -6 < z < 6 for a grid N r ,N$,N z = 100x6x200. 
Results of simulations are shown in Fig. 11 The 
color background shows pressure distribution (left 
panel) and density (right panel). The lines show 
corresponding analytical solutions. 


5.2. Test of the diffusivity module 



Figure 11: Test of the ideal MHD module of the 3D version 
of code with the Low-type solution test at different times. Top 
Panels: Pressure (left panel) and density (right panel) distribu¬ 
tions (color background) and analytical Low-type solution (solid 
lines) with grid resolution N r ,N^,N z = 100 x 6 x 200 at scaling 
parameter a = 2.15, Bottom Panels: scaling parameter a = 2.4. 



Here we test the diffusivity module of the code. 
For this we switch off the hydrodynamic fluxes and 
the corresponding right-hand side terms in Eqs. 


and 37 and integrate these equations numerically 
in the region 0 < r < r\, -r\ < z < n with the 
rectangular region 0 < r < r 0 , -r 0 < z < r 0 excised 


Figure 12: The top panels show contours of the azimuthal com¬ 
ponent of the magnetic field B magnetic flux function T and 
the azimuthal velocity from left to right. The solid lines in¬ 
dicate the numerical solutions while the dashed lines show the 
analytical solutions given by Eqs. |52||54| The bottom panels 
show the absolute errors in the calculation of the corresponding 
functions. 
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Figure 13: An example of hydro simulations of the interaction of the counter-rotating disks. Top Panels show the case where the inner 
and outer parts of the disk are rotating in the opposite directions. Bottom Panels show the case where a clump of gas rotates in the 
direction opposite to the rotation of the disk. The color background shows the density distribution, vectors show matter flux, and 
the white bold line separates the regions of opposite angular velocities which are marked with the “+” and signs. Simulations are 
performed in 2d cylindrical coordinates with the grid N r x N, = 150 x 250. From|Dyda et al. 1 12015). 


from simulations. 
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d'V 
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d 

~dt ~dr 
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1 3 doj 1 

d 

i 

r 2 dr 
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, + dz 

r* 


(48) 


(49) 


= 0 . 


(50) 


The solutions for the magnetic flux function 'F 
and azimuthal velocity are symmetric about the 
equatorial plane while solutions for the azimuthal 
component of the magnetic field B<p are antisym¬ 
metric. All solutions go to zero at the symmetry 
axis 


'F| r= o = 0, Bp | r= o - 0, v^| r =o - 0. 


We can find particular solutions of Eqs. 48 50 by 
the method of variable separation in spherical coor¬ 


dinates. For t] = ], Eq. 48 in spherical coordinates 


is 


d'V d lx ¥ dll dT) 
dt dR 2 d9 \ sin 6 d9 ) 


(51) 


For testing we choose the following solutions of Eq. 

|5ll 


^ = C sin 2 6 ^^ _ Ssimg j^ + 10(f 0 - f) j 

CsinOcosOt 6(to -1)\ 

B * = — r —l 1 + 

Csin0/ 5sin 2 0\/ lO(to~t)\ 

v * = ~i^V - —A 1 + — 


(52) 

(53) 

(54) 


Grids 

maxIAT 

max|AB^| 

max Av'^l 

50 x 100 

3.68 x 10 ~ 2 

3.19 x ur 2 

7.59 x 10 -2 

100 x 200 

1.01 x 10~ 2 

6.51 x 10~ 3 

1.60 x 10~ 2 


Table 1: The maximum absolute error in calculations of the 
T, and v$ obtained in the diffusivity tests. 


where C and to are constants, R = Vr 2 + z 2 is the 
spherical radius, and 9 is the polar angle. 

The equations for the magnetic flux function and 
the ^-components of the velocity and magnetic field 
are integrated numerically using our 2.5D code 
(with the switched-off hydrodynamic fluxes and 
zero right hand side terms). We use the same grid 
geometry as in the main simulations; that is, cylin¬ 
drical coordinates with a non-equidistant grid in 
both directions. 

For the test, we take a grid of N r ,N z = 50 x 100, 
and also a finer grid with N r ,N z - 100x200. We also 
take the inner and outer radii (in the /--direction) 
of the simulation to be ro = 1 and /y = 15.6. For 
these boundaries we set V F, B^ and to be equal 
to the analytical values determined by Eqs. 52|5T 


The initial conditions correspond to the exact solu¬ 
tions of Eqs. |52||54 at t — 0. The constants C and to 
in these equations are chosen so that the solutions 
at the inner boundary are of the order of unity. We 
use C = 0.01, to = 50 and integrate the equations 
up to t\ = 2to = 100 . 

The top panels of Fig. 12 show the result of the 
diffusivity test for the 0 -component of the magnetic 
field, the vector potential, and the 0 -component 
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of the velocity. The solid lines show the numer¬ 
ical solution of Eqs. 48|50 


which are obtained 


by integrating the equations using our code un¬ 
til t - 100. The dashed lines show the exact so- 
52||54 


lutions given by Eqs. 


Evidently, the nu¬ 


merical and analytical solutions are nearly identi¬ 
cal. The bottom panel shows the absolute errors 
Ah' = T'num - M J ex;ict for each function. One can see 
that the error is very small everywhere, taking into 
account the fact that the maximum value of all the 
functions, V P, v# » 1. The error increases to¬ 
wards the star because the value of the functions 
also strongly increases towards the star. We per¬ 
form similar simulations and analysis for finer grid 
resolutions 100 x 200. Table [l] shows the maximum 
absolute errors, max|'F num -'F exact |, max|B™”'-B" flC '|, 
and max|v™ m - v e ™ c, \. One can see that for the finer 
grid (dimensions 100 x 200), this error is 2.1 and 
3.9 times smaller compared to the coarsest grid (di¬ 
mensions 50 x 100). This confirms the convergence 
of the numerical solution towards the analytical so¬ 
lution. Similar comparisons for the B (b component 
give factors 2.2 and 4.1, which also show conver¬ 
gence. Similar comparisons for the v (! , component 
give factors 2.2 and 4.1, which also show conver¬ 
gence. 

6. Examples of Astrophysical Simulations 

The codes discussed in this paper have already 
been used to study a variety of different astrophys¬ 
ical problems. Below, we show a few examples of 
calculations with the code in both 2D and 3D ge¬ 
ometries. 

6.1. Counter-rotating accretion disks: 2.5D hydro 
simulations in cylindrical coordinates 

The hydrodynamic version of the code in 2.5D 
cylindrical coordinates was used to study counter¬ 
rotating accretion disks (Dyda et al.[ 2015). The 


simulations show that the interaction between two 
counter-rotating regions of the disk strongly de¬ 
pends on the value of the viscosity coefficient a 
which varies from a = 0.02 up to a = 0.4. The 
simulations show that in counter-rotating disks, 
the accretion rate may be increased by a factor of 
~ 10 2 - 10 4 . 

The problem is set up such that the inner region 
of the disk orbits the star in one direction while the 
outer part orbits in the opposing direction (see Fig. 
13 top panels). The simulations were performed 


for a viscosity parameter of a = 0.1. Initially, the 
boundary between the counter-rotating regions of 
the disk is sharp. When the simulations start, the 
oppositely directed regions of the disk start inter¬ 
acting due to viscosity, leading to annihilation of 
the angular momentum at the interface. This mat¬ 
ter, which has lost some of its centrifugal support, 
begins to accrete rapidly due to the unbalanced 
gravitational force. This leads to a hundred-fold 
increase of the accretion rate compared with a disk 
which does not have a counter-rotating component. 
The simulations also show that the rcj> component 
of the stress tensor, T rih dominates the angular mo¬ 
mentum transport in the disk. 

We also study the case where a clump of counter¬ 
rotating matter is allowed to interact with the top 


of the disk (see Fig. 13 bottom panels). We find 


that a part of the disk matter loses angular momen¬ 
tum through viscous interactions with the counter¬ 
rotating clump and is quickly accreted onto the star. 
Analysis shows that in this case, three components 
of the viscosity tensor r^, t :i!i and r r - contribute to 
the angular momentum transport. 

6.2. MRl-driven accretion onto a magnetized star: 
2.5D MHD simulations in cylindrical coordi¬ 
nates 

The code has also been utilized to study accre¬ 
tion onto a magnetized star from a turbulent ac¬ 


cretion disk in 2.5D cylindrical coordinates (Ro- 
manova et al.[ 2011). A sample of the results are 


shown in the right panel of Figure [14] The turbu¬ 
lence present in the disk is driven by the magneto- 
rotational (MRI) instability (e.g., Balbus & Hawley 
1991). Simulations were performed for different 


magnetospheric radii r m and different initial values 
of the plasma parameter p = %np!B 2 in the disk, 
which varied from p = 10 to p = 10,000. The dy¬ 
namical range of initial densities between the disk 
and corona was 10 4 . Inside the magnetosphere, the 
dipole magnetic field varies as ~ r~ 3 and the mag¬ 
netic pressure varies as ~ r~ 6 . In the regions of high 
magnetic pressure, the matter pressure and density 
become very low, thereby limiting the simulation 
timestep and slowing the calculation. To overcome 
this difficulty, we require: p > /3 floor x B 2 /8n, where 
Pfloor = 10 3 ; in other words, we require that the 
matter pressure does not fall below a certain frac¬ 
tion of the magnetic pressure. In the grid cells 
where the matter pressure becomes smaller than 
this value, it is adjusted to be p = /? floor x B 2 /%n. We 
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Figure 14: From left to right: (1) The cylindrical grid N r x N z = 270 x 432 cells used in simulations of the magnetospheric accretion, 
compressed near the axis and in the disk; (2) The same grid, rarefied by a factor of 10, for demonstrative purposes; (3) An example 
of a simulation of the disk-magnetosphere interaction, where color shows the density distribution, and black lines are magnetic field 
lines; (4-6, top panels) An example of the disk-magnetosphere interaction in the case of an MRI-driven turbulent disk, the bottom 
panel shows the accretion rate. From Romanova et al. ( 20111. 


then recalculate the density p = ( p/s) l/y to be ther¬ 
modynamically consistent. In cases of large mag¬ 
netospheres, this leads to slightly enhanced mat¬ 
ter density and pressure near the star. However, 
this does not change the physics of the magneto¬ 
sphere because the magnetic pressure is still 10 3 
times higher than the matter pressure, and there¬ 
fore the physics of the disk-magnetosphere interac¬ 
tion is described correctly in the sense that the mag¬ 
netosphere is strongly magnetically-dominated. In 
other parts of simulation region, no floor density or 
pressure is added. 

The grid is compressed towards the disk mid¬ 
plane and towards the z-axis with the goal of hav¬ 
ing higher grid resolution in the disk and near the 


star. The left panel of Figure 14 shows a typical 


grid used in this study with N r = 270 cells in the 
radial direction and N z = 432 cells in the axial di¬ 
rection. The number of grid cells covering the disk 
in the vertical direction (in the middle of the disk) 
is about 200, and the MRI-unstable modes are well- 
resolved even in cases of a small-scale turbulence. 

For this study, the viscosity or diffusivity terms 
are not enabled. The disk matter is loaded onto the 
magnetic field lines of the magnetosphere due to 
reconnection or due to the small numerical diffu¬ 
sivity which is present in the code. 

6.3. Propeller outflows from the disk-magnetosphere 
boundary: 2.5D MHD simulations 
Another problem which has been studied using 
this code is that of disk accretion onto a magne¬ 


tized star rotating in the propeller regime (Lii et 


ah] 2014). In the propeller regime of accretion, 
the magnetosphere of the star rotates more rapidly 
than the inner disk; the angular momentum im¬ 
parted on the disk by the magnetosphere can push 


the disk away or drive outflows from the disk- 


magnetosphere boundary (e.g., Illarionov & Sun- 


yaev 1975; Lovelace et al. 1999). As in the previous 


study (Sec. 6.21, the angular momentum transport 
in the disk is mediated by the MRI-driven turbu¬ 
lence. 


The large centrifugal barrier at the disk- 
magnetosphere boundary inhibits direct accretion 
onto the star. Rather, the disk matter accumulates 
at the disk-magnetosphere boundary until it com¬ 
presses the magnetosphere inward of the corota¬ 
tion radius and accretion can proceed. The result 
is a burst of accretion in which part of the mat¬ 
ter accretes onto the star and part is launched as 
a short-lived magneto-centrifugal wind. There is 
also a secondary low-density wind which flows out 
along field lines anchored to the surface of the star. 
In the latter component of the winds, the matter is 
launched rapidly, resulting in low density regions 
near the star. 


In order to prevent the density (and simula¬ 
tion timestep) from becoming too small, we im¬ 
plement a floor density condition. That is, we fix 
the minimum density in the simulation region to 
be pfi 00 r = 0.01 Pcorona = 10 6 , effectively adding a 
small amount of matter whenever the density drops 
below a certain threshold. We track the amount 
of matter that is added, and find that it is much 
smaller than the matter flux in the disk and in the 
wind. We also implement the matter pressure floor 
condition for matter pressure inside the magneto¬ 
sphere (see Sec. 6.2). 


As a word of caution: for relatively large val¬ 
ues of /?A oor , the matter which is added to main¬ 
tain the floor pressure can have a tendency to flow 
out of the super-Keplerian magnetosphere due to 
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Figure 15: An example of a 2D MHD simulation of propeller driven outflows performed in cylindrical coordinates with non-uniform 
grid of dimension N r x N z = 250 x 432 cells. Panels a and b: energy flux density carried by the magnetic fields, F gf, and by matter, F g m . 
The lines show the magnetic field lines. Panel c: the contours show the poloidal matter flux, pv p , and the red arrows show the poloidal 
velocity vectors. Panels d and e: angular momentum flux density carried by the magnetic fields, F/,/ z , and matter, F/, m -. FromlLii et al.l 
l2M4l l. 


centrifugal forces, thereby forming an artificial out¬ 
flow. However, for low enough values of /3f\ OOI (up 
to jSfioor = 3 x 10 5 ), no artificial outflows from the 
inner magnetosphere were observed. 


Fig. 15 shows matter flux and also energy and 
angular momentum carried by the matter and the 
magnetic field (see details and description of fluxes 
in |Lii et al.| 2014). The axisymmetric grid is in 
cylindrical (r, z) coordinates with mesh compres¬ 
sion towards the disk and towards the z-axis such 
that there are a larger number of cells in the disk 
plane and near the star. In the models presented 
here, we use a non-uniform grid with dimension 
250 x 432 cells, corresponding to a grid that is 66 
by 140 stellar radii in size. At r = 20, the num¬ 
ber of grid cells which cover the disk in the vertical 
direction is about 200. 

At the disk-magnetosphere boundary, 3D insta¬ 
bilities may cause the disk matter to penetrate 
through the magnetic field lines into the magne¬ 
tosphere. Since these instabilites are intrinsically 
three dimensional, they are not modeled in the 
current axisymmetric 2.5D simulations. To mimic 
these processes, we implement diffusivity near the 
star using the diffusivity module described in Sec. 


4.2.2 To prevent modifications to the overall disk 
accretion rate, the diffusivity is only enabled near 
the disk-magnetosphere boundary We study the 
disk-magnetosphere interaction in this region for 
different diffusivity parameters ad = 0.01,0.03,0.1. 
We find that for larger values of ad, the disk mat¬ 
ter penetrates through the magnetosphere more 
rapidly causing most of the accreting matter to be 
ejected into an outflow 


6.4. Migration of a planet in the magnetized disk: 

2D MHD simulations in polar coordinates 

A version of the code in 2D polar coordinates 
has been used to model the migration of a planet 
through a magnetized disk (Comins et al. 2015, 
in prep.). In a test version of this problem, we 


followed the approach by Fromang et al. (2005) 
which used two different codes to solve this prob¬ 
lem (FARGO and another code). Our initial con¬ 


ditions are almost identical to those taken by Fro¬ 


mang et al. (2005): the planet’s mass is 5 Earth 


masses; the disk has a constant surface density dis¬ 
tribution; the initial magnetic field has only an az¬ 
imuthal component, I hi,, which scales with radius 
as Bq, = B,f,o(ro/r) k where k = 0,1,2; and lastly, the 
disk is locally-isothermal. We take y = 1.01 for 
the adiabatic index to mimic a locally-isothermal 
disk. The size of simulation region is R min = 0.4 and 
Rmax = 4.9, and the grid N r = 480 and = 1200. 
The grid is stretched in the radial direction such 
that the grid cells are almost rectangular in all parts 
of the simulation region. The magnetic resonances 
are resolved at 12 grid cells in the radial direction. 

Our simulations confirm the results obtained by 
Fromang et al. ( 2005] ). Namely, that if the plasma 
parameter [i - %np/B 2 is of order unity, then the 
magnetic resonances can be more important than 
the Lindblad resonances (which lead to inward mi¬ 
gration) and the migration can be stopped if the 
radial gradient of the azimuthal magnetic field is 


sufficiently large. Fig. 16 shows the density dis¬ 


tribution where the magnetic resonances are ob¬ 
served as azimuthally-stretched density depletions 
near the planet (see left panel). The middle panel 
shows that in the region of magnetic resonances, 


21 


























































Figure 16: Modeling a planet in a magnetized disk in 2d polar coordinates. The planet excites slow magnetosonic waves at the magnetic 
resonances which are seen as low-density rings in the left panel and an enhanced azimuthal field in the middle panel. The right panel 
shows the linear distribution of the surface density and magnetic field in the region of the resonances. Simulations performed with the 
grid N r = 480, and A(j = 1200. From Comins et al. (in prep.). 


the magnetic field is enhanced. The right panel 
shows a radial cut of the density and magnetic 
field along the planet’s radial vector. We find that 
at k = 0 the migration is inward; in contrast, for 
steeper density profiles (k = 1 and k = 2) the planet 
migrates outward. This is in accord with the re¬ 
sults obtained by Fromang et al. (2005|, who find 
that the inward migration consistently slows and 
reverses direction as k increases. 


6.5. Migration of a planet in non-magnetized lam¬ 
inar disk: 3D hydro simulations in cylindrical 
coordinates 

A version of the 3D code in cylindrical coordi¬ 
nates has been used to model the migration of a 
planet in a 3D hydrodynamic disk. The cylindri¬ 
cal simulation region extends between R min = 0.4 
and R max = 4.9 and Z max = +0.5. The number of 
grid cells in each direction is N r = 232, N$ = 480, 
and N z = 80. The simulations show that the planet 
excites two waves in the disk at the Lindblad res¬ 
onances and migrates inward due to angular mo¬ 
mentum transfer from the planet to the disk. The 
density perturbation is largest in the equatorial 
plane, and falls off away from the equatorial plane. 


Fig. 17 shows an example of simulations which 


show typical waves in the disk. 


choose the adiabatic form of the energy equations. 
The codes pass standard tests and were used for 
solving different astrophysical problems. 

The simulations show that our entropy- 
conserving HLLD solver performs satisfactorily 
in many astrophysical applications, such as studies 
of low-density, high-velocity outflows in propeller- 
driven winds. The ideal MHD module, viscosity 
module and diffusivity module were tested in 
different astrophysical situations. 
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